Synchronization in A Carpet of Hydrodynamically Coupled Rotors with Random 

Intrinsic Frequency 
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We investigate synchronization caused by long-range hydrodynamic interaction in a two- 
dimensional, substrated array of rotors with random intrinsic frequencies. The rotor mimics a 
flagellated bacterium that is attached to the substrate ("bacterial carpet") and exerts an active 
force on the fluid. Transition from coherent to incoherent regimes is studied numerically, and the 
results are compared to a mean-field theory. We show that quite a narrow distribution of the intrin- 
sic frequency is required to achieve collective motion in realistic cases. The transition is gradual, 
and the critical behavior is qualitatively different from that of the conventional globally coupled 
oscillators. The model not only serves as a novel example of non-locally coupled oscillators, but 
also provides insights into the role of intrinsic heterogeneities in living and artificial microfluidic 
actuators. 
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Introduction. - Collective oscillations of active ele- 
ments are observed in a variety of physical, chemical, and 
biological systems far from equilibrium. Numerous stud- 
ies have been devoted to the mutual entrainment of os- 
cillators that have different intrinsic frequencies 0, 0] ■ A 
class of phase oscillators with global (or mean-field) cou- 
pling have enjoyed deep theoretical understanding [lj-Q) 
while a myriad of unresolved problems still remain on the 
behaviors of locally Q and non-locally 0, 0] coupled os- 
cillators. In particular, knowledge about synchronization 
caused by long-range interactions is quite limited 0-Q> 
although they are ubiquitous in Nature in the form of, 
e.g., gravitational, electromagnetic, elastic, and hydro- 
dynamic forces. 

Biologically important examples of long-ranged syn- 
chronization are provided by swimming microorganisms 
that are interacting hydrodynamically, such as sperm 
flagella beating in harmony (lfj| - |T3 |. and metachronal 
waves in cilia 15-l{j. Both of these oscillatory elements, 
flagella and cilia, are driven by molecular motors embed- 
ded in the cell surface, and interact through the viscous 
environment (water). In order to describe their synchro- 
nized motion, several theoretical models have been pro- 
posed USE SEMI. While all of these initial studies 
are focused on a homogeneous set of identical systems, it 
will be important to consider the role of disorder, as real 
biological motors possess intrinsic heterogeneities that 
could affect the collective dynamics. 

An example of collective yet heterogeneous dynamics is 
found in a bacterial carpet, which is recently introduced 
as a new type of microfluidic device [20j . The assem- 
bly is composed of a dense monolayer of bacteria that 
are attached to a solid substrate by their bodies (heads) . 
Their flagella (tails) , on the other hand, can freely rotate 
in the fluid and are orientationally ordered by hydrody- 
namic interaction, to generate coordinated fluid motion. 



Evolution of correlated regions is observed, but the order- 
ing remains partial. Observations of irregular and slowly 
varying flow structures ('whirlpools' and 'rivers') sug- 
gest the presence of heterogeneity in the configurations 
of the rotors 21] . Fabrication of more efficient microflu- 



idic pumps could be achieved through understanding and 
controlling the heterogeneity. 

Recently, we have proposed a generic model of hydro- 
dynamically coupled rotors arrayed on a substrate, and 
studied the collective dynamics of uniform elements |22j . 
In this letter, using a variant of the model, we address 
the effect of random intrinsic frequencies on synchroniza- 
tion. To be concrete, we consider a simple and idealized 
model of bacterial carpets, where a flagellated bacterium 
is mounted to each rotor nearly (but not completely) ra- 
dially. Depending on the mounting angle of the flagellum, 
it generates a torque that drives the rotor in either clock- 
wise or counterclockwise direction. We assume that the 
mounting angle and hence the intrinsic frequency of the 
rotor are randomly distributed. By varying the degree 
of randomness, we study the transition from coherent to 
incoherent regimes numerically. Our results suggest that 
synchronization of the rotors, and hence collective pump- 
ing of the fluid, requires a quite narrow distribution of the 
mounting angle in realistic cases. In order to understand 
the transition behavior, we apply the mean-field theory, 
which is originally developed for global coupling, to our 
long-ranged system. While we obtain a fair agreement 
between theory and simulation for the synchronization 
threshold, the transition is shown to be more gradual 
than in globally coupled oscillators. 

Model. - We consider an array of rotors positioned 
on a square lattice of grid size d. Each rotor has a thin, 
freely-rotatable arm on the tip of which a flagellated bac- 
terium is mounted. The bacterium consists of a spherical 
bead of radius a (body) and a thin tail that lies horizon- 
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FIG. 1: Schematic picture of the rotor (top view). The active 
force Fi exerted by the i-th rotor on the fluid is deviated by 
a fixed angle Si from the radial direction. The reaction force 
drives the rotor at the intrinsic frequency w, = F sin 5i/(b, 
where b is the radius of rotation and £ the viscous drag coef- 
ficient of the bead. 



tally (flagellum). Motion of the bead is constrained on 
a circular orbit of radius b located at height h from the 
substrate, which we take to be the a;y-plane. The posi- 
tion of the i-th bead is thus given by = r oi + he z + bm 
where roi is its base position on the square lattice and 
rii = (cos fa, sin^, 0) is the unit vector that gives the ori- 
entation of the arm via its phase fa = fa(t). The velocity 
of the bead reads Vj = fati where ij = (— sin fa , cos fa , 0) 
is the unit vector tangential to the trajectory. We as- 
sume that the active force fj exerted by the rotor on 
the fluid has a constant magnitude F, and makes a fixed 
angle Si (measured clockwise) from the radial direction; 
Fi = F(cos5iTii — sin5,tj). See Fig. [T]for the configura- 
tion. The reaction force —Fi on the rotor arm gives the 
driving torque = Fb sin Si and the intrinsic frequency 
uii = F sin Si I Q>, where C = dnrrja is the viscous drag co- 
efficient. The mounting angles Si's are assumed to have 
the Gaussian distribution 
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with the standard deviation So- 

We assume that the rotors are widely spaced so that 
a,b,h -C d. Then the velocity field of the fluid created 
by the active forces is given by v(r) = J^i G*(r — ■ Fi, 
where G(r) = (3h 2 /2nr]) ■ rj_rj_/|r| 5 , rj_ = (x, y, 0) is 



the asymptotic expression of the Oseen-Blake tensor [23 1 
for h/d <C 1. The rotor's angular velocity is given by 
uji + v(rj) • ti/b, or, more explicitly, 



dfa 
~~dt 



too sin Si 
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Here, w = F/Qb, ry = Vi-Vj = |ry|(cos%,sin%), and 
7 = Qv 1 /nd 3 = 6irah 2 /d 3 is the dimensionless coupling 
constant. For a real bacterial carpet, a ~ h ~ lpm and 
d ~ 10/im give the rough estimate 7 ~ 10 -2 . 

Numerical Simulation. - We implemented the model 
on a L x L square lattice and numerically integrated 



eq. ([2]) by the Euler method. We assumed the peri- 
odic boundary condition and computed the velocity field 
every time step in the Fourier space. We set 7 = 0.1 and 
varied the angle deviation <5 as the control parameter. 
The system size used was L = 128 for most of the results 
shown below, while L = 32, 64 and 256 are also used to 
check finite-size effect. Starting from random initial con- 
figurations of tf>i(t = 0), the system reached a dynamical 
steady state by the time t = 1 x 10 4 /wo- The statisti- 
cal data shown below are taken from the time window 
1 x 10 4 < u t < 2.5 x 10 5 . 

We plot the orientational order parameter S = | (n) \ = 
as a function of Sq in Fig. O Also shown is the 
standard deviation (STD) of the actual frequency f2^ = 
{4>i) normalized by the STD of the intrinsic frequency 
u>i, Q = y/(n 2 )/(u 2 ). Note that S = 1 and Q = in 
the fully synchronized state and 5 = and Q = 1 in 
the desynchronizcd limit. As we increase So, S and Q 
slowly converge to the desynchronized limit. While the 
change in the orientational order parameter is sharper for 
a larger system size, the frequency deviation has little in- 
dependence. For comparison, we also show the results 
of the mean-field theory, which will be explained in the 
next section. 

In Fig. [31 we plot the distribution function of the ac- 
tual frequency f2 normalized by the STD of intrinsic fre- 
quency, for different values of 8q. The distribution con- 
sists of a sharp delta-function like peak at ft — and 
broad symmetric tails for > and < 0. For So < 3°, 
most of the rotors are coherent and contribute to the 
center-peak. For Sq = 10°, the distribution is close to 
that of the intrinsic frequency, while the center peak still 
remains, the above data suggest that the synchroniza- 
tion transition in this system is more gradual than that 
found in a globally coupled system, and it is difficult to 
locate the transition point exactly. 

On the other hand, the variance of the order parameter 
Var(S) = (S(t) 2 ) - (S(t)) 2 as a function of S (Fig. Ha)) 
has a peak near Sq = 6°, suggesting that there is a subtle 
balance between synchronization and desynchronization. 
We will call this the threshold angle and denote by 5th- 

In FigUfb), we plot the temporal correlation function 
of the order parameter C s {t) = (S(t+t')S(t'))~(S) 2 . We 
find an oscillatory behavior with long correlation time at 
So = 6°. The presence of a characteristic period is also di- 
rectly observed in the plot of S(t) in Fig. 0|c). Although 
the origin of the oscillation is beyond the scope of the 
present paper, a preliminary study shows that the oscil- 
latory behavior at the threshold angle is a unique feature 
resulting from the long-ranged nature of the interaction. 

We also plot the orientational correlation function 



G n (|r|) = ([n(r + r')-n]>(r')-n]} 



(3) 



in Fig[5] Here the angular brackets mean taking average 
over r' as well as the azimuthal angle of r. For 4 < So < 
6°, we observe an exponential decay of the correlation 
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FIG. 2: Transition behavior for system size L — 32, 64, 128, 
and 256 with comparison to the mean-field theory (MFT). (a) 
Orientational order parameter S = |(n}| = K^)! versus So- 
(b) Normalized frequency deviation Q — \J (Q 2 )/{lj 2 ) versus 
So versus So- 
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FIG. 3: Distribution of the normalized actual frequency 
Q/y/Jufi) as functions of So. For comparison, the distribu- 
tion of the intrinsic frequency Wj is also shown. 
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FIG. 4: (a) Variance of the order parameter Var(S') as a 
function of So- Fluctuation is most enhanced at So — 6° . (b) 
Auto-correlation function of the order parameter Cs(t). An 
oscillatory behavior is prominent at t5o = 6°. (c) Time series 
S(t). 



over a wide distance. For Sq > 6°, on the other hand, the 
correlation is short-ranged and decays more slowly than 
exponential. The qualitative change in the correlation 
function gives another support of the above estimate of 
the threshold angle, S t h = 6°. 

Mean-field theory. - Now we apply a mean-field ap- 
proximation to our model. When the coupling is weak 
(7 <C 1) and the force angle is small (<5o <C 1), we can 
neglect the Sj 's in the RHS of eq. because they give 
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FIG. 5: Orientational correlation function G n (r). Nearly 
exponential decay of correlation is observed for So < 6°. 



O(7<$o) contributions. We also replace the interaction 
kernel oc rr/r 5 by its angular average, as a result of which 
the cosine term in the RHS of eq. @ is dropped. This 
gives the phase equation in the familiar form, 



dt 



Tj) sin( 



f) 



(4) 



with u>i — uj sinSi and G(r) = Sjuiod 3 /4nr 3 . The dis- 
tribution of the intrinsic frequenc y is given by P u (u>) = 
\dS/duj\P(S) = P(sin _1 {lo /ujq)) / \J<jJq — u 2 . Now we ap- 
ply the mean-field ansatz which was originally proposed 
by Kuramoto [H for global coupling: 



Re v 



(5) 



Here, the amplitude R and the phase 6 of the order pa- 
rameter are assumed to be constant in space and time, 
which would be justified if the interaction is sufficiently 
long-ranged. The isotropy of G(r) allows us to assume 
6 = without loss of generality. Using this we can rewrite 
eq. (1) as 



dt 



Rs'mc 



(6) 



Equation (3) allows a stationary (ci = 0) solution if and 
only if \tJi\ < R. The rotors satisfying this condition have 
the actual frequency O = and are called the coherent 
group. The phase of a rotor belonging to this group is 
given by 



3111 



(7) 



where the principal value of the inverse-sine function 
should be chosen so that \<p\ < w/2. The phase distri- 
bution n(<j)) of the coherent group reads 



n{<f>) = P u (u) 



doj 



= P w (R sin <j>) ■ R cos < 



(8) 



The rotors with |w(r)| > R, on the other hand, form 
the incoherent group. The actual frequency of a rotor 
belonging to this group is given by 



2tt 



2n A A dt 



Jo d< t> 



-R 2 . 



(9) 



The phase distribution n'(<f>) = n'(<p; tOj) of an incoherent 
rotor is proportional to the frequency it comes to </>: 



i'(<p;uJi) = Clcir 1 = C\uj t -Psinc 



(10) 



with the normalization factor C = y/wf — R 2 /2ir. 

Now we replace the factor e l ^ j in the RHS of eq. §5§ 
by its ensemble average as 



R = Y J G(v-v'){e^). 



(11) 



The average is the sum of the contributions 
from the coherent and incoherent groups, 
(e**)coh = Zi) 2 d(t>n{(j>)e^ and (e^) mcoh = 
f , >R du)P u (u)) J_ n dc/m'(4>; w)e l< *'. The latter van- 
ishes because of the symmetry of P w (u;), and the former 
with eq. flS) yields 

R = RGqJ(R), (12) 

,tt/2 

J(R) = / dcjiP^R sin c/>) cos 2 <j>, (13) 

J -it/2 

with Go = J2 r G( r )- This is a self-consistent equation 
for the mean-field amplitude. Expanding the integral as 
J(R) = (tt/2) [P w (0) + Pj(0)P 2 /8 + 0(P 4 )] , we obtain 
the critical coupling strength 



Goc 



7rP u (0) 



(14) 



for the synchronization transition (at which a non- 
vanishing solution R appears). For a square lattice, 
we have Go = 9.03 • 37Wo/4-7r. Also we have Pj(0) = 
1/ (v2ftOJo$o) ■ Putting these together into eq. (fl"4"|) . we 
obtain the critical angle 



S 0c = 1.357. 



(15) 



In the simulation we used 7 = 0.1, which gives 5q c = 
0.135 (rad) = 7.73 (deg). This value is not very far 
from the numerically obtained threshold angle 5th = 6 
(deg). Both the simulation and theory suggest that a 
quite narrow distribution of the mounting angle is re- 
quired to achieve coordinated motion in a real bacte- 
rial carpet. On the other hand, the mean-field theory 
predicts a sharp transition, in contrast to the gradual 
crossover observed in the simulation. Near So Cj the orien- 
tational order parameter decays as S oc v?0c ~ <^0i while 
the normalized frequency deviation linearly approaches 



5 



to the desynchronized limit as 1 — Q oc Sq c — 5q. The dis- 
tribution of the actual frequency Pq(TI) for S < 5q c has 
three distinct peaks, one at uj = (the coherent group) 
and two symmetric peaks for u) > and u) < (the in- 
coherent group). The central peak vanishes for 5 > So c - 
These behaviors are qualitatively different from the nu- 
merical results shown in Figs. 2 and 3. In order to explain 
the unconventional transition behavior, it is necessary to 
develop a theoretical framework that incorporates spa- 
tial fluctuations, which is an interesting problem for the 
future. 

Conclusion. - The synchronization transition caused 
by long-range hydrodynamic coupling is shown to be 
more gradual than for the global coupling. The threshold 
angle for the crossover is estimated, and is not far from 
the mean-field estimate for the transition point. It sug- 
gests that a very narrow angle distribution is required to 
achieve correlated motion in a real bacterial carpet. We 
believe that this work sheds some light on the delicate is- 
sues involved in hydrodynamic synchronization and hope 
that it stimulates further theoretical and experimental 
studies of the rich behavior of such systems. 
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